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O ' ABSTRACT 

<N ; 

Observations have shown that the Eddington ratios (the ratio of the bolometric luminosity 
to the Eddington luminosity) in QSOs/active galactic nuclei (AGNs) cover a wide range. In this 

■ paper we connect the demography of AGNs obtained by the Sloan Digital Sky Survey with the 
accretion physics around massive black holes and propose that the diversity in the Eddington 
ratios is a natural result of the long-term evolution of accretion disks in AGNs. The observed 

• accretion rate distribution of AGNs (with host galaxy velocity dispersion a ~ 70 — 200 km s -1 ) 

^\ , in the nearby universe [z < 0.3) is consistent with the predictions of simple theoretical models 

in which the accretion rates evolve in a self-similar way. We also discuss the implications of the 
results for the issues related to self-gravitating disks, coevolution of galaxies and QSOs/AGNs, 
I/"") , and the unification picture of AGNs. 

■ Subject headings: black hole physics — accretion, accretion disks — galaxies: active — galaxies: 
evolution 

o ■ 

CZJ . , . 

■ 1< Introduction 

• i— i , The past decade has seen great progress on measuring the mass of massive black holes (BHs) in both 

nearby normal galaxies and AGNs/QSOs (e.g., Magorrian et al. 1998; Richstone et al. 1998; Wandel et al. 
1999; Kaspi et al. 2000). With reliably measured BH masses in AGNs/QSOs, their Eddington ratios (i.e., 
the ratio of the bolometric luminosity to the Eddington luminosity) can be estimated. Compared to the 
bolometric luminosity (or the accretion rate, given the mass-to-energy conversion efficiency e), the Eddington 
ratio (or dimensionless accretion rate) is a normalized parameter and more related to the accretion physics. 
Observations have shown that the Eddington ratios in AGNs/QSOs cover a wide range (at least 2-3 orders of 
magnitude) with an upper limit of about 1 (Woo & Urry 2002; McLure & Dunlop 2004). It is natural to ask 
what determines the diversity of the Eddington ratios. The simplest attempt to answer this question would 
be by first asking whether the diversity of the Eddington ratios is a consequence of evolution of accretion 
processes onto BHs. 

The questions raised above on the origin of the diversity of the Eddington ratios in AGNs/QSOs may 
have important implications for understanding the evolution of the AGN/QSO luminosity function (e.g., 
Schmidt & Green 1983; Boyle et al. 2000). Generally the number density of AGNs/QSOs at a given red- 
shift and luminosity/dimensionless accretion rate should be jointly determined by both the number den- 
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sities of galaxies whose nuclear activities are triggered at earlier time and the follow-up nuclear luminos- 
ity/accretion rate evolution of individual AGNs/QSOs. The history of the triggering rate is related to the 
formation and evolution of galaxies (e.g., major mergers of galaxies) and may have a strong cosmic evolution 
(e.g., Kauffmann & Haehnelt 2000). The follow-up nuclear luminosity /accretion rate evolution of individual 
AGNs/QSOs is directly related to the evolution of the accretion disk around the BH and the fueling process 
onto the disk in their host galaxies. In this paper we extract the evolution of the accretion processes from 
the observational distribution of accretion rates in AGNs. This provides constraints on theoretical models 
of accretion disk evolution and also illuminates the origin of the diversity of the Eddington ratios. 

To extract the information on the nuclear luminosity/ accretion rate evolution from observations, sta- 
tistical methods involving a large sample of AGNs are required since (i) a single AGN may only represent 
one specific period in a prolonged phase of nuclear activity. A large sample of AGNs with different ages will 
span all phases of this activity and allow us to extract information about evolution; and (ii) in addition to 
age, other physical parameters may be important in determining how AGNs evolve and a statistical method 
may help to clarify these. In § 2, we describe our method of extracting the evolution of accretion processes 
phcnomcnologically from the observational distribution of the accretion rates in nearby AGNs (see also the 
related methods in Yu & Lu 2004a, b). We argue that results for nearby AGNs are not strongly affected 
by the uncertainty in the nuclear activity triggering history of these objects, because the evolution of the 
triggering rate at low redshifts {z < 0.5) is probably rather weak. Recent progress by the Sloan Digital Sky 
Survey (SDSS) on the demography of AGNs (http://www.mpa-garching.mpg.de/SDSS/; Kauffmann et al. 
2003b; Heckman et al. 2004) makes the application of the method practical. In § 3, we apply the method 
to a large sample of SDSS AGNs and get the phenomenological constraints on the evolution of accretion 
processes. We find that the constraints are consistent with theoretical expectations of the accretion disk 
evolution models shown in § 4. This consistency suggests that the diversity of Eddington ratios in AGNs is 
probably a natural result of the evolution of the accretion disks around massive BHs. In § 5, we also discuss 
the implications of the results for the long-term evolution of accretion disks, AGN unification models, and 
coevolution models of galaxies and QSOs/AGNs etc. Our conclusions are summarized in § 6. 



2. Methodology 

In this section we analyze how the accretion rate evolution of an AGN is incorporated into the accretion 
rate distribution of AGNs in the local universe. 

Observations have revealed that the central BH mass in nearby normal galaxies is tightly correlated with 
the velocity dispersion of the bulge components of host galaxies (Ferrarese & Merritt 2000; Gebhardt et al. 
2000; Trcmaine et al. 2002). In this paper we take this correlation as a boundary condition for mass growth 
of BHs in AGNs, that is, the final BH masses of nearby AGNs after the nuclear activity, Mf, is assumed 
to follow the same correlation with velocity dispersion a, since mass growth of BHs comes mainly from the 
nuclear activity phase (e.g., Yu & Tremaine 2002; Fabian 2004; Marconi et al. 2004). Here the velocity 
dispersion a (or the gravitational potential field in galactic scales) is assumed to not change significantly 
during and after the nuclear activity phase. This assumption has also been implicitly used in many physical 
models to interpret the correlation between BH mass and velocity dispersion in nearby normal galaxies 
(e.g., Silk & Rees 1998; Blandford 1999; Fabian 1999; Adams et al. 2001; Burkcrt & Silk 2001; King 2003; 
Murray et al. 2005; see also more discussions on the BH mass versus velocity dispersion relation in AGNs 
in § 5 and in Yu & Lu 2004b). We describe the correlation through a probability distribution function of 
logM/ at a given er, P(\ogMf\a), and P(log Mf\a)dlog Mf gives the probability that the logarithm of the 
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BH mass of an AGN with velocity dispersion a is in the range log Mf — ► log Mf + d log Mf . The mean 
log M f at a given a is denoted by (logM f \a)[= 8.13 + 4.021og(cr/200kms- 1 ); Tremaine et al. 2002]. The 
intrinsic scatter of P(logM/|<r) around (logM/|cr) is small enough that it is difficult to distinguish from 
measurement errors, but is less than 0.27 dex (Tremaine et al. 2002). Given the bolometric luminosity 
of an AGN Lbob the mass accretion rate is given by Mboi = L^ \/(ec 2 ), where c is the speed of light. 
In practice, the [OIII] A5007 line luminosity Ljoiii] is often used as a tracer of AGN activity (e.g., see 
Hcckman et al. 2004); and the bolometric luminosity of an AGN may be estimated from its observed [OIII] 
line luminosity, given its bolometric correction /boi = £boi/£[onrj- We may also define the accretion rate 
through the [OIII] line luminosity by Mom] = /b o i^[0in]/( ec2 )7 where / bol is a factor representing the 
average bolometric correction between -Lpm] and Lboi and given through log/£ ol = (log/bol)- The average 
of the logarithm of the bolometric correction (log/ b oi) may be obtained by giving the probability distribution 
function P(logL[ ni] l^boi), which is defined so that P (logiroiii]|iboi)^[oiiil 1S the probability that the 
logarithm of the [OIII] line luminosity is in the range log L [oni] ~^ l°g^[oni] + ^l°g^[oiii] at a given Lboi- 
In general, Mboi and Mjohj] are not the same since the bolometric corrections are not the same for all the 
AGNs. 

Given an AGN with nuclear activity triggered at cosmic time ti and with final BH mass Mf, we describe 
the evolution of its accretion rate by Mbo\{T, ti, Mf), where r = t — ti is the physical time passed since the 
triggering of the nuclear activity. The total time that the AGN will spend with the logarithm of its accretion 
rate in the range log Mboi —* logMboi + dlogMboi at cosmic time ti = t(z = z{) < t < t(z = 0) = t is given 
by (see eq. 15 in Yu & Lu 2004a): 

T(M hoh U,Mf)dlogM hoi 

^ (Mboi In 10)d log Mboi ^ 

where (k = 1,2,...) are the solutions of the equation M ho \{Mf,T,ti) = M bo i that satisfy t\ < tj + 
T k < t , and T(M hol ,ti, Mf) = if no such solutions exist. In this paper the rcdshift z ; is set to 
0.3 and AGNs with z < zi arc taken as nearby AGNs. We define the function (p(M^ i, Mf, a, t) so 
that (f>(Mb \, M f , a,t)dlog M^oidlog M fdloga gives the comoving number density of AGNs at redshift z(t) 
with the logarithms of the accretion rate, the final BH mass, and the velocity dispersion in the ranges 
log Mboi — * log Mboi + d log Mboi, log Mf — ► log Mf + dlogMf, and logc — > loga + dlogcr. Given the evolu- 
tion of the nuclear activity triggering rate N(U, Mf, a), which is defined so that N(U, Mf, a)dtid\og Mfdloga 
is the comoving number density of BHs with host galaxy nuclear activity triggered at cosmic time U — ► U+dU 
and with the logarithms of its final BH mass and velocity dispersion in the ranges log Mf — > logMj+dlogM/ 
and log a — > loga + dlogo, we have the time integral of (f>(Mb \,Mf,a,t) given by (cf., eq. 10 in Yu & Lu 
2004a) 

fto 

<j>(M hol ,Mf,a,t)dt 



I 

Jti 



= f N(ti,Mf,a)T(M hoh ti,Mf)dti 
Jo 

= N(M hoh Mf,a)T(M hoh Mf), (2) 

where 

_ Jo" N(U,M f , a)T(M hoh U,M f )dt. 
Jl T{M hoU U,Mf)dt 



N(M hoh M f ,a) ee Jo K ;;;zl ( 3 ) 
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and 



T(M hoh M f ) = I ° T(M hoh ti,M f )dU. (4) 



By multiplying equation (2) by P(log L[om] |iboi) and then integrating it over log Mboi and logM/, we have 



\ 



where 



$("l[OIII],<7, t)dt 

= j dlogM f J dlogM hol P(logL [OUI] \L hol )N(M hoh M f ,a)T(M hoh M f ), (5) 

$(ra [om] ,c7,i) = J d\ogM f J d log M bol P (log L [OIIl] \ L bol ) (f>(M hol ,M f , a, t), (6) 

m [0 ni] = M[oiii]/M E dd[M}(c7)] is the normalized accretion rate, M' f (a) = io< lo « A/ / l CT > , M E dd[M f (cj)} = 
LEdd[M'j:(a)]/ (ec 2 ), and L^dd[Mf(a)] is the Eddington luminosity of a BH with mass Mf (a). Note that here 
the definition of the normalized accretion rate is a little different from the conventionally defined Eddington 
ratio in that -Mjom] may not be the same as M bo i, and the Eddington luminosity used here is for the 
average final BH mass, rather than for the BH mass that is powering the nuclear activity in an AGN. If 
the bolometric correction of each AGN in the sample and its BH mass that is powering the nuclear activity 
are obtained from future observations, we may construct the methods for the distribution of the Eddington 
ratios as above and perform a consistency check on the results. 

We can further define the time-averaged accretion rate and velocity distribution function of nearby 
AGNs (with z < zi) by 

,. . tf* <S>(m [olllh <j,t)dt 

®(rn [OIU] ,a) = -I — . (7) 

Below we assume that the triggering rate of nearby AGNs N(ti, Mf,a) = n(ti, a) P (log Mf\a), where 
n(ti,a) — J dlogMf N(ti, Mf,cr) is independent of Mf. By further assuming that n(ti,a) is independent of 
ti, we have N(Mboi, Mf,a) = N(ti,Mf, a); and then by combining equations (5) and (7), we have 

^(m [0 ui],<y) 

= n(a) J dlogM f P(logMf\a) j d log M hol P(log L [om] \L hol ) T( ^'^ /) ■ (8) 

Below we also assume that A^boifV, ti, Mf) is independent of the cosmic time ti, and according to equations 
(1) and (4), we have 

T(M M ,M,) = W'"~" . M "'° 10 *, (9) 

,„ Jmin{0,ti-T k \ dMbo\(T,Mf)/dT\ 



where Tk (k = 1,2, ...) are the solutions of the equation A4-t, Q i(Mf, t) = M bo i that satisfy < < to. If 
all the solutions < ti (k — 1,2, which is satisfied by the assumed models and its best-fit parameters 
shown in Fig. 2 below), we have T(M bo b-Mf) oc t — ti, ^(n^oni], cr) — $(m[oin] , a , t) at z(t) < z\, and 
<J>(to[oiii], ct, t) is independent of t. In reality AGNs may have a high triggering rate at high redshifts, 
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however, AGNs triggered at high redshifts may have become faint enough at the current epoch to move out 
of the range of the accretion rates of interest in this paper. For example, using the accretion rate evolution 
model assumed below (eqs. 10 and 11) and the best-fit parameters obtained in § 3, only those AGNs whose 
nuclear activities are triggered at z < 0.5 may have an accretion rate high enough at z < z\ to be within 
the range shown in Figure la below. As an additional check, we have assumed that the triggering rate is 
proportional to the major merger rate of galaxies, e.g., n(t i} a) oc [1 + z{ti)] a with a = 0.51 ± 0.28 (Lin et 
al. 2004), and our results are not significantly affected. 

It is plausible to assume that the growth of a BH involves two phases after the nuclear activity is 
triggered on (Yu & Lu 2004a; see also Small & Blandford 1992). In the first (or "demand limited") phase, 
there is plenty of material to fuel BH growth. However, not all of the available material can contribute to the 
BH growth immediately and the BH growth is limited by the Eddington luminosity. We assume that the BH 
accretes with the Eddington luminosity for a period tj and denote the BH mass at cosmic time t j = U + tj 
by Mj. We assume the mass-to-energy conversion efficiency e to be a constant. Thus, the accretion rate in 
the first phase increases with time as follows: 

-Mboi(r) = M Edd (Mj) exp {^^j > < t < tj, (10) 

where ts p = 4.5 x 10 7 e/[0.1(l — e)]yr is the Salpeter time [the time for a BH radiating at the Eddington 
luminosity to e-fold in mass; we set e/(l — e) = 0.1 below]. With the decline of material supply, the BH 
growth enters into the second (or "supply limited" ) phase and the nuclear luminosity is expected to decline 
below the Eddington luminosity. Below we assume that the accretion rate in the second phase evolves as 
follows: 

M bol (r) = M Edd (Mj) ( T ~ ^ + TJ ) 7 , t>tj, (11) 

where -Mboi( T ) K at t — tj 3> t d and t d is the characteristic transition timescale from the first phase 
to the power-law declining stage of the accretion rate. We assume that the two phases appear only once for 
an AGN. The final BH mass obtained after the nuclear activity is given by 

/•OO 

M f = Mj + / (1 - e)M ho i(T)dt 

Mj, (12) 



1 + 



TD 



(7- l>Sp_ 

where the efficiency e is assumed to be the same as the efficiency in the first phase. Here we do not consider 
that the accretion process may transit to a radiativcly inefficient mode when the accretion rate becomes low 
enough (i.e., when the Eddington ratio is < 0.001; Narayan & Yi 1994; Blandford & Bcgelman 1999). Such 
rates are in any case below the accretion rate range studied in this paper (see Fig. 1). For simplicity, below 
we assume that the timescales tj and td are independent of Mf. 

According to the accretion rate evolution model assumed above and equations (l)-(6), ^(rhroni], cr) is 
roughly constant with logmjom] if most AGNs arc in the first phase and is proportional to rh~^l if most 
AGNs arc at the late power-law declining stage of the accretion rate [here P(log Mf \a) and P(log £[om] l-^boi) 
are assumed to follow Gaussian distributions; see also § 3]. If the scatter of the distributions P(logM/|cr) is 
sufficiently small and the bolometric correction /boi is roughly a constant, m[oni] has a maximum given by 

m max = [l + tTD/TBpVfr-l)]- 1 (13) 

when the AGN is transiting from the first phase to the second phase. 
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3. Application to observations 

The SDSS has obtained a large sample of 33,589 narrow-line AGNs in the local universe (http://www.mpa- 
garching.mpg.de/SDSS/Data/agncatalogue.html). Details of the sample selection and the derivation of the 
AGN properties used in this paper [Xromii ct, [i* (the stellar surface mass density)] are given in Kauffmann 
ct al. (2003a, b) and Heckman et al. (2004). In this section we fit the accretion rate evolution models assumed 
in § 2 (eqs. 10 and 11) to these observed AGNs. (Note that Seyfert I AGNs are not included in this sample 
because for these objects the AGN itself may have a significant impact at least on the estimate of the velocity 
dispersion. Possible effects on our conclusions if including Seyfert I galaxies will be discussed at the end of 
this section.) 

We consider AGNs with L[oiii] > £[om],min = 10 6 L Q , a > cr m i n = 70kms _1 , and (i, > 3 x 
10 8 M Q kpc~ 2 . A sample of fainter AGNs is incomplete because at larger distances, the fiber samples more 
of the stellar light of the host galaxy and weak emission lines become difficult to detect. The lower limit 
on velocity dispersion is set by the instrumental resolution of the SDSS spectrograph. The galaxies with 
lower stellar surface mass densities are mainly disk-dominated (Kauffmann et al. 2003a), and their velocity 
dispersions may be substantially affected by the disk components. (Note that the quantity £[om] used here 
is not extinction-corrected. The extinction effect may be implicitly included in the scatter of the bolometric 
correction used in the calculations below. See more discussions in Heckman ct al. 2004). We get the obser- 
vational distribution $o' bs (TOrom], ct) (j = 1, 2, ..; k = 1, 2...) by counting the AGNs in such bins: (i) the bin 
of the velocity dispersions is [log <Tj , log Oj+i] , where log(<7j/er m i n ) = (j — 1) A log a and the bin size A log a = 
0.05dcx. (ii) The bin of the normalized [OIII] line accretion rate is [logm[oni],fc,logm[om],fc+i], where 
lo g(^[oni],fcM[oiii],mm) = (k- l)Alogm [0 m], logm [0 iii],mm = log L [OIII] , min - log L Edd ((M f )) + log/ bol , 
and the bin size AlogTO[ OIII ] = Alog a x d log [Mj(a)]/d log a ~ 0.20dex. Each observed AGN i has a 
weighted factor l/V maXt i, where y max ,i is the volume over which the object would have been detectable 
(Schmidt 1968). The quantity ^^(frironi], a ) is obtained by summing up the weighted factors (^2 i l/V maXt i) 
of AGNs in each bin, and its error <5$ obs is estimated by (^ 1/Knax i) 1 ^ 2 ( see points and their dotted error 
bars in Fig. 1). Note that the realistic error may not be well represented by the estimate above if the number 
of AGNs in a bin is small. Below in the fit we do not use the points with only one (open triangles in Fig. 1) 
or two (open squares) AGNs in their bins. In each velocity dispersion bin, the falloff of $o' bs (?7i[oiii] , a) at 
the low-accretion-rate end (see crosses in Fig. 1) is caused by the luminosity cut at £[oni] = ^[oni],min and 
these crosses are not used in the fit, either. 

We estimate the model distribution of ^mod^iom]) a ) by urst using equation (8) to obtain $mod("i[om]> ct), 
and then binning the distribution and normalizing it by setting ^ fe ^^(wfoiii]) a ) = J2k ^obU^IOiii]) CT )' 
where the sums are taken over the bins with $o' bs (?ri[oni]i ct) ^ (but excluding the bins of crosses and open 
points shown in Fig. 1; similarly for the sum in \ 2 below). We assume that P(logMf\a) follows a Gaussian 
distribution with an intrinsic scatter of 0.2 dex around (logM/|cr) (the value of the scatter is set with some 
degree of arbitrariness, since as mentioned in § 2 the intrinsic scatter is difficult to distinguish from mea- 
surement errors but is less than 0.27 dex; Tremaine et al. 2002. See also discussions below). The distribution 
of P(log L[om] |iboi) is assumed to be independent of Lboi and also follow a Gaussian distribution with 
(log(l//boi)) = (log(L [0 iii]/£boi)) ~ log(l/3500) and a scatter of - 0.38 dex (Heckman et al. 2004). We use 
the x 2 statistic to fit the model distribution to the observational distribution of ^'^(^[oiiili u ) (j = 1; 2, ..; 
k = 1, 2, ...). The best-fit parameters are obtained by minimizing \ 2 = Sj.fc[(log $^ d — log 3>^ s )/i51og ^o'^J 2 - 
We use P x 2 to denote the probability of x 2 being higher by chance if the data were drawn from the model. 
Our a priori definition of an acceptable fit is one with P x 2 > 0.01. For different t//ts p = 5,1,0.1,0, we 
obtain the best fits for parameters (7,td/ts p ). We find that if the observational distribution only with 
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a\ < a < trio — 200 km s -1 (the range of the black hole mass M'^{a) over the nine bins of the velocity 
dispersion is 2.0 x 10 6 -1.3 x 1O 8 M0; see Fig. la) is included in the fit, all the best fits are acceptable with 
P x 2 ~ 0.02; if the observational distribution only with o\ < a < erg ~ 176 km s -1 is included in the fit, 
the best-fit parameters are not affected significantly but with an improved P x 2 ~ 0.07 — 0.10; and if adding 
the observational distribution with higher velocity dispersions a > cio (see Fig. lb) in the fit, no fits are 
acceptable (since the distribution may be complicated by some other factors, not solely by the accretion 
evolution, which will be further discussed below). 

For the observational distribution with <J\ < a < cio, the contour plots of confidence levels for the 
best-fit parameters are shown in Figure 2. The obtained best-fit parameters (see crosses in Fig. 2) arc 
insensitive to tj and the average of the best fits over four cases of t//ts p is (7,td/ts p ) ~ (1.26,3.1). We 
take the errors of the fit parameters as ~ ±0.1 and 5(td/ts p ) — ±1, which roughly correspond to the 68 
per cent confidence contours (i.e., x 2 — Xmin — 2.30) shown in Figure 2. We show the model distribution 
*mtd( T "[oiii]:°') U = 1,2,. ..9) obtained with (t//t Sp , 7, t d /t Sp ) ~ (1,1.26,3.1) as solid lines in Figure la. 

As shown in Figure la, the observational accretion rate distribution at 10~ 3 < m[ ni] < 0.1 are well 
consistent with a power-law distribution, and the parameter 7 is mainly determined by the slope of the 
distribution (— 1/7). If assuming an exponentially declining form of the accretion rate evolution in the 
second phase of the accretion model instead of the power-law declining form, all fits will be unacceptable. 
This can be easily understood, since the exponentially declining form predicts a flat accretion rate distribution 
(see eq. 1), as illustrated by the horizontal dashed line in Figure 1 (assuming no significant cosmic evolution 
of the nuclear activity triggering rate). 

The best-fit parameter tjj is mainly determined by the observed turnover of the power-law distribution 
of the accretion rates at the high-rate end (tojoiii] ~ 0.1; see also m max in eq. 13). In Figure 1, although 
there exist AGNs with logmpin] > logm max ~ —1.1, the existence of these AGNs does not necessarily 
mean that they are accreting at a super-Eddington rate with Mbo\ > A^Edd(A^BH) [here MEdd(^fBH) is the 
conventionally defined Eddington accretion rate for the BH mass Mbh that is powering the nuclear activity 
in the AGN, not the final mass Mf after the nuclear activity]. In the accretion evolution model used here, 
we always have Mboi < M^^Mer) ■ The existence of m[ OIII ] > m max is caused by the non-zero scatter of 
the Mf — a relation and the bolometric correction. The best-fit parameter td may be affected by the values 
of the scatter in the Mf — a relation and the bolometric correction. For example, if the intrinsic scatter 
of P(\ogMf\cr) is 0.27 dex, the best fits are (7,T£>/ts p ) ~ (1.24,3.5); if the intrinsic scatter is zero, no fits 
are acceptable, but the best fits for the observational distribution with o\ < a < 09 are acceptable with 
(7, t d/t~s p ) — (1.23, 2.1). Precise understanding of the evolution of accretion processes would require precise 
measurements of this scatter (see also Yu & Lu 2004a). 

As mentioned above, the observational sample used here does not include Seyfert I galaxies. The fraction 
of type II AGNs in the SDSS is a weak function of AGN luminosity, decreasing from ~ 60% to ~ 30% over 
the range Lpm] ^ 10 5 5 - 10 8 5 L© (Hao ct al. 2005; Hcckman et al. 2004). Thus, the fraction of AGNs with 
high accretion rates (e.g., at tojoiii] ~ 10 _1 — 1; see Fig. 2) could be underestimated here, which will mainly 
affect the estimate of td- After including type I galaxies, the best-fit parameter td would decrease and 7 
would increase slightly, resulting in an increase in m max (eq. 13). If we assume that the dependence of the 
fraction of type II AGNs on accretion rate is the same as the dependence on luminosity, after including type 
I AGNs, we find that the distributions (at the high accretion rates) shown in Figure 1 move upward by at 
most log 2 ~ 0.3 dex. The value of td/ts p would then have to change from 3.1 to 2.2 in order to compensate. 

The best-fit td above is obtained by using the number distribution of AGNs as a function of accretion 



- 8 - 



rate. Generally, the timcscale related to the length of the nuclear activity may also be constrained by 
comparing the number density ratio of AGNs to their remnants (e.g., Yu & Tremaine 2002). The predicted 
ratio of the number density of nearby AGNs to nearby inactive galaxies is affected by uncertainties in the 
triggering history of AGNs at high redshifts and hence does not help us obtain a more precise constraint. 

As mentioned above, if adding the observational distribution with high velocity dispersion a > ctio — 
200 km s™ 1 (see Fig. lb) in the fit, no fits are acceptable. A full explanation for their detailed distribution 
would need to involve at least the following four factors, (a) These objects generally have relatively low 
normalized accretion rates, and thus their accretion rate distribution could be affected by the change of 
the efficiency e due to the transition from the thin-disk accretion mode to the radiative inefficient accretion 
mode, (b) Host halos or galaxies of high-mass BHs generally have different major merger histories than those 
of low-mass BHs. The AGNs with high velocity dispersions generally have high BH masses, low normalized 
accretion rates, and probably long ages of the nuclear activities, and thus their distribution could be more 
likely to be affected by the cosmic evolution of the nuclear activity triggering history, (c) The low accretion 
rate in these objects (e.g., m < 10 -3 ) may be simply caused by Bondi accretion. Note that the observational 
accretion rate distribution at the low accretion rate end shown in Figure lb has indicated a steeper increase 
than the slope — 1/7. (d) The assumption that td are independent of Mf may be violated for a wide range 
of Mf. 



4. Comparison with accretion disk evolution models 



The process that removes the angular momentum of gas in a galaxy and drives it into the very central 
region is still not well understood. It may require a hierarchy of mechanisms, handing off to one another as 
the material moves inward. Here we consider the "last stop" of the "handing-off" . We hypothesize that some 
mechanism has moved gas to the central parsec in the galaxy, where it will then diffuse to the central BH 
through viscous transport (or specifically by the magneto-rotational instability; e.g., Balbus & Hawley 1998) 
in an accretion disk. The surface mass density, S, of an accretion disk with a central, dominant gravitating 
point mass and with viscosity v evolves as a function of time r and radius R through the equation 



or. 



3__a_ 

RdR 



R^±(^R^) 



(14) 



If the viscosity is of the form v(R, S) oc R m T, n (where m and n are constants), this equation has solutions 
of self-similar form (Pringle 1981) as follows: 



S(A,T) 

So 



R_ 

Ro 



(15) 



where the values of 77 and £ are determined by the exponents m and n in the function v{R, £) and also by 
the detailed boundary conditions of the accretion disk. So, -Ro, and ro are scale variables with arbitrary 
dimensions and they satisfy To = Rq/v(Rq, So)- The self-similar evolution of accretion disks has also been 
calculated using numerical simulations (e.g., Lin & Pringle 1987; Cannizzo, Lee & Goodman 1990; Pringle 
1991). For an accretion disk with negligible torque at both the origin R = and its outer boundary (so 
that its total angular momentum is conserved during its evolution), the rate at which the mass of the disk 
decreases is given by a power law of the form M d (T) = ^(J 2ttT,RcLR) oc t 32+170+21, j where a and b are 
the exponents in the opacity law n(p,T) — n p a T~ b (e.g., Cannizzo, Lee & Goodman 1990). If a = b = 
(Thomson opacity), we have r\ = -15/16, £ = -3/8, f(x) = (28)~ 3/2 x- 3 / 5 (l - x 7 / 5 ) 3 / 2 in equation (15), 
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Fig. 1. — Distribution of the normalized accretion rate mjoni] m nearby AGNs (z < 0.3). The difference 
between mpni] and the conventionally defined Eddington ratio is described in § 2. Different colors represent 
different velocity dispersion bins [log <Tj, log<7j + i] (J — 1,2,...) with o\ = 70kms _1 and bin size Aloger = 
0.05 dex. In panel (a), from top to bottom, the colors black, red, green, blue, cyan, magenta, black, blue, and 
green correspond to j — I, 2,. ..,9, respectively (cio — 200 kms -1 ); and in panel (b), the colors red, cyan, 
magenta, black, blue, and green correspond to j = 10, II, ...15, respectively. The sum Sj in the normalization 
of the distribution is over j = 1,2, ..9 in panel (a) and over j = 1,2, ...15 in panel (b). The observational 
distribution obtained from nearby Seyfert II galaxies are drawn as points (triangles, crosses, or squares) at 
the middle point of each log(TO[om]) bin and their error bars are shown by dotted lines. For open triangles, 
only one AGN is in each bin; and for open squares, only two AGNs are in each bin. The falloff of crosses at 
the low-mpiii] end is caused by the luminosity cut at £[oni] = £[om],min- The binned model distributions 
(only for j — 1,2,. .9 in panel a), obtained with (t//ts p , 7, td/ts p ) — (1,1.26,3.1), are connected by solid 
curves (for view clarity, the binned distributions are not drawn in the form of a histogram). The observational 
distribution at the accretion rates 10 -3 < TOpm] < 0.1 is well consistent with a power-law distribution. The 
tilted dashed lines are the reference lines for the power-law distribution with a slope of — 1/7. The horizontal 
dashed lines are the reference lines illustrating a flat distribution predicted by an exponential declining form 
of the accretion rate evolution in the second phase of the accretion model. See details in § 3. 



- 10 - 



1.4 - 



1.2 - 



1.4 - 



1.2 - 



1 1 1 1 1 II 1 1 1 1 1 1 1 1 1 1 

; T iAsp= 5 ; 


i i i 1 1 1 1 1 i i i i 1 1 1 1 1 
" T lA S p =1 

/ 


: T ./%=o^ ; 

1 1 


1 1 1 1 III 1 1 1 1 1 1 1 1 1 

: T ./ T s P =o^ ; 

1 1 



- 1.4 



- 1.2 



- 1.4 



- 1.2 



1 



10 



Sp 



1 

T D A 



10 



Sp 



Fig. 2. — Contours of the confidence levels for the best-fit parameters 7 and td- The velocity dispersions of 
the AGNs used in the fit are in the range [a\, <7io] — [70, 200] kins" 1 . The best-fit parameters are shown by 
crosses. The contours enclose the 68.3%, 95.4%, and 99.7% confidence regions on them in each panel. See 
details in § 3. 
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and the mass accretion rate evolves as 



|Ai(T)|= (WWry^ 



70 \T 0/ 

where Md,o = (28)~ 3 / 2 ^-RqSo, and the diffusion timescale of the disk 



ro = K^k) = ((U - L6) x 10 yr Vo.3 - ipc; Vio 7 m g 



f- 

Vo.i 



1 /• M dfi N 2/3 
10 7 M R 



(17) 



where a is the commonly used parameter related to the viscosity in the standard accretion disk model, and 
i?o may be taken as the size of the accretion disk or the location of its reservoir. The value of To is sensitive to 
i?o and the reference value of Rq in equation (17) is roughly consistent with the observations of the accretion 
disks in several nearby AGNs (e.g., NGC4258, Miyoshi et al. 1995; NGC1068, Lodato & Bertin 2003; and 
NGC 3079, Kondratko et al. 2005). If a = 1 and b = 7/2 (Kramers opacity), we have \M d \ oc t~ 5 / 4 . 

Note that the best-fit parameter 7 = 1.26±0.1 for AGNs with a c± 70 — 200 kms -1 is consistent with the 
exponent (-1.18 or -1.25) of r in M d expected from the self-similar solutions above. The diffusion timescale 
To in equation (17) is also roughly consistent with td obtained from observations. These consistencies 
between the observations and the theoretical expectations above are an encouraging result. They support 
the hypothesis that the observed diversity in the Eddington ratios of AGNs results from an evolutionary 
sequence of accretion processes. They also suggest that the accretion disks of most Seyfert galaxies (with 
c c± 70 — 200 km s -1 and m > 10~ 3 ) arc currently in the self-similar phase or are transiting toward this phase. 
It is interesting to note that the self-similar evolution of disk evolution is also supported by observations of 
accretion disks around T Tauri stars (Hartmann et al. 1998). 



5. Discussions 

Despite the encouraging consistency above, the self-similar solution of the accretion disk evolution shown 
in § 4 is based on highly simplified accretion disk models. For example, the effects on the disk evolution 
by disk winds or by infalling material that is continuously deposited onto the disk are not considered. The 
role of possible disk instabilities is not considered, either. More complete investigations on the long-term 
evolution of accretion disks are needed. 

There are also further issues to be considered: 



• Self- gravitating disks and star formation: With the best-fit parameters obtained for the accretion 
evolution in § 3, the disk is self-gravitating in its early evolutionary stage and the central BH mass 
becomes dominant in its late evolutionary stage. We note that this evolutionary pattern of the disk 
mass is probably compatible with current observations of the accretion disks in a few AGNs. For 
example, NGC3079 (Kondratko et al. 2005) and NGC1068 (Lodato & Bertin 2003) have Eddington 
ratios close to 1, which may correspond to the first "demand-limited" phase or a transition to the 
self-similar evolution of the second phase. Their accretion disks appear to be thick with masses larger 
than or comparable to their central BH masses. NGC4258 has an Eddington ratio < 0.01 (Gammie et 
al. 1999 and references therein), so it may be in the late stage of self-similar evolution, and its accretion 
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disk is thin with a mass that is much smaller than its BH mass (Miyoshi et al. 1995). The evolution of 
a self-gravitating disk may not exactly follow the simple prescription adopted in this paper. Recently, 
Fromang et al. (2004) simulated the evolution of self-gravitating tori around a massive BH and found 
that they quickly evolve to a dual structure with an inner thin Keplerian disk fed by an outer thicker 
self-gravitating disk (sec also Bcrtin & Lodato 1999; Lodato & Rice 2004, 2005). Further investigation 
of the long term evolution of self-gravitating disks is crucial for understanding the transition to the 
self-similar phase. In addition, in a self-gravitating disk, star formation may be important because of 
gravitational instabilities that develop in the disk (Goodman 2003; Levin 2003). However, the growth 
of the central BH will consume a substantial fraction of the disk mass. The fraction of the disk mass 
that goes into stars rather than the BH is worthy of further investigation. 

• The Mf — a relation in normal galaxies and the Mbh — c relation in AGNs: As mentioned in § 2, 
in our model the Mf — a relation in normal galaxies is used as a boundary condition for BH mass 
growth in AGNs, and BH masses in AGNs are lower than the expectation from the Mf — a relation. 
In observations some AGNs with sufficiently low Eddington ratios may be included in normal galaxies. 
As a self-consistency check for our model, here we show that the best-fit parameters for BH growth 
obtained in § 3 do not contradict the tightness of the Mf — a relation. For example, using equations 
(11) and (12), for (7,td) = (1.26,3.1), when the Eddington ratio of an AGN declines to below 10~ 3 , 
the logarithm of the ratio of its central BH mass to its final mass is — 0.18dex < log(MBH/M/) < 0, 
which does not exceed the scatter of the Mf — a relation used in the model. In addition, considering 
that most of the supplying material may be blown away as outflows when a BH accretes material via 
the ADIOS (adiabatic inflow-outflow solutions) mode at low accretion rates (Blandford & Begelman 
1999), the difference between the BH masses in AGNs with low Eddington ratios and their final masses 
will be smaller. 

For the best fits (7, t d ) = (1.26, 3.1), BH masses in most AGNs (a ~ 70 - 200 kms" 1 ) with Eddington 
ratio in the range 1-0.1 are about 0.08-0.3 times their final mass Mf. (Note that here at a given 
velocity dispersion and luminosity, most AGNs are in the second evolutionary stage, rather than in the 
first stage in which the nuclear luminosity exponentially increases.) However, by studying 15 dwarf 
Seyfert I galaxies (a ~ 30 - 80 kms" 1 ), Barth et al. (2005) found that the BH masses in these AGNs, 
Mbh (estimated from the empirical relations among the BH mass, the nonstellar 5100A continuum 
luminosity, the broad- line region radius, and the broad- line width in Kaspi et al. 2000), arc consistent 
with the Mf—cr relation and these AGNs have Eddington ratios in the range 0.1-1. Currently it is hard 
to answer whether the evolution model in our paper contradicts the observations or not, because the 
answer should take into account at least the following uncertainties, (i) Usually, BH masses Mbh used 
to derive the empirical relations in AGNs are estimated through the reverberation mapping technique 
with some assumptions on the kinematics of the broad-line emission clouds (e.g., Wandel et al. 1999; 
Kaspi et al. 2000). The reverberation mapping technique could either systematically overestimate or 
underestimate the BH mass by a factor of 3 or so (Krolik 2001). (ii) It is possible that the sample in 
Barth et al. (2005) biases toward detection of AGNs with the highest mass black holes for a given a. 
(iii) The constraint on td obtained in § 3 may not be strict (e.g., because of the exclusion of Seyfert 
I galaxies). If td is smaller (e.g., ~ lrs p ), the BH masses in most AGNs with Eddington ratios in the 
range 1-0.1 will be closer to their final mass (e.g., ~ 0.2M/ — 0.4M/). (iv) The dwarf galaxies in Barth 
et al. (2005) are in a velocity dispersion range different from the SDSS AGNs and may not follow the 
same BH growth parameters. 

In addition, we note that the BH masses in AGNs measured through the reverberation mapping 
technique by Onken et al. (2004) are calibrated to follow the Mf — a relation in normal galaxies by 
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multiplying them by a factor of 1.8. Onkcn ct al. (2004) argue that the calibration factor could be 
caused by a nonisotropic velocity distribution of the broad-line emission clouds; however an alternative 
explanation could be that there exists an offset between the Mbh — o relation in AGNs and the Mf — a 
relation in normal galaxies, which just represents BH mass growth during the nuclear activity phases 
(Yu & Lu 2004b). It is worthwhile to study whether such an offset really exists, which would provide 
strong constraints on BH growth models and the physical origin of the Mf — a relation. 

• Binary BHs (BBHs) and coevolution of galaxies and QSOs/AGNs: In the current cocvolution model 
of galaxies and QSOs (Kauffmann & Haehnelt 2000), the nuclear activity of galaxies are triggered by 
major mergers of galaxies. Mergers of galaxies may form a BBH if there exists a massive BH in each of 
the merging galaxies (Begelman, Blandford & Rees 1980; Yu 2002; Volonteri, Haardt & Madau 2003). 
The evolution of the accretion disk around a massive BBH may be different from that around one BH 
(e.g., because of different boundary conditions of the disk) and depend on the orientation of the BBH 
orbital plane and BH spins. For example, studies of the time-dependent evolution of the accretion 
disk around the primary of a binary star with the outer disk radius tidally truncated by the secondary 
also yield self-similar solutions but with a steeper time dependence of the accretion rate Md oc r -2 ' 5 
(Thomson opacity) or r~ 3 - 3 (Kramers opacity) (Lipunova & Shakura 2000; see also a different self- 
similar solution for an external disk around a binary star by Pringle 1991). The consistency of the 
slope of the accretion rate function with the expectation of our simple model may suggest that galaxy 
mergers have not been sufficiently efficient in faint AGNs to form a large population of BBHs, or that 
galaxy triaxiality or gas dynamics has caused BBHs to merge very quickly. 

The AGN luminosity function, defined by the comoving number density of AGNs per unit luminosity 
at a given redshift, is proportional to J m _1 4>(m, o~)d\oga with m cx L / L^^M'^a)]. The luminosity 
function at the faint end, if fitted by a power-law form cx L@ , is expected to have the slope < — 1 — 1/7 
(i.e., < -1.85,-1.80,-1.40,-1.30 for 7 = 1.18,1.25,2.5,3.3, respectively), where the symbol < is 
expected instead of the symbol ~ because the break of the power-law accretion rate distributions 
(m ~ m max ) at the high-rate end corresponds to different break luminosities for AGNs with different 
velocity dispersions, and thus the faint end of the luminosity function consists of not only AGNs 
with accretion rates lower than the break rate but also some with rates higher than the break rate, 
which may steepen the slope. According to the fit for the accretion rate distribution of nearby AGNs, 
the expectation for their luminosity function at the faint end, < —1.8, is consistent with the fit 
for the observational luminosity function by Hao et al. (2005). Observations of the QSO luminosity 
function indicate a relatively flat slope ~ —1.45 ± 0.03 at the faint end (Richards et al. 2005), which 
might suggest the existence of BBHs in a significant fraction of QSOs since —1.8 < < —1.3. It is 
worthwhile incorporating the evolution of BBHs and the evolution of accretion disks in models for the 
coevolution of galaxies and QSOs/AGNs. One could then make joint predictions for the evolution of 
the QSO/ AGN luminosity function as well as their accretion rate distribution function for comparison 
with observations. 

• Possible difference between AGNs and bright QSOs: The luminosity evolution of bright QSOs (Lbol > 
10 46 ergs _1 ) is constrained from observations through a relation between the BH mass function in 
nearby normal galaxies and the time integral of the QSO luminosity function in Yu & Lu (2004a). The 
characteristic timescale for the luminosity of the QSO to decline during its second evolution phase was 
shown to be very short (< 0.3rs p ; Yu & Lu 2004a), significantly shorter than the timescale obtained 
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for AGNs in this paper. 1 This difference between AGNs and bright QSOs is slightly smaller if we 
include type I AGNs, but it still persists (see discussion in § 3). In Yu & Lu (2004a), the luminosity in 
the second phase is assumed to decrease exponentially and not according to a power law. If the QSO 
luminosity evolution is assumed to have the same evolutionary form as described in equations (10) and 
(11), an even smaller td (and/or larger 7) would be required. According to equation (17), To (or td) 
is highly sensitive to the location of reservoir of accreting materials (oc R ) but depends weakly on 
Mbh and M4. A smaller td for QSOs implies that the disk size or the location of the material that 
fuels the BH in QSOs is not larger than in AGNs; this is surprising because the BHs in QSOs are one 
or two orders of magnitude larger than the BHs of the AGNs that we have studied in this paper. This 
may provide interesting constraints on the fueling mechanisms that remove the angular momentum of 
materials at galactic scales and cause them to circularize at the parsec scale. A larger 7 might suggest 
the existence of BBHs in QSOs (see also the discussion above) . 

• Unified Models of AGNs: adding the effect of an evolving accretion disk In current unified models 
of AGNs, orientation effects are proposed to explain the diversity of AGNs/QSOs (Antonucci 1993). 
The accretion rate and the Eddington ratio have also been proposed as important parameters for 
understanding why AGNs are so diverse (Boroson 2002). In this paper, we connect the diversity of the 
Eddington ratio in AGNs with the evolution of accretion disks, and we suggest that different kinds of 
AGNs can be partly unified by considering the time sequence of the disk evolution. 

— Lack of a torus in very weak AGNs: The current AGN unification model requires the existence of 
a torus in an AGN, but observations show that it does not exist in M87 (Whysong & Antonucci 
2004) . Here we propose that the lack of a torus may be explained as an evolutionary consequence 
of the accretion process. At an early stage, the outer disk may be thick and clumpy and play the 
role of a torus. At very late stage, most of the outer disk (or the torus) is very likely to have 
diffused away and the emission from the rest of the torus may be too weak to detect. M87 has a 
low Eddington ratio < 10~ 3 (Di Matteo et al. 2003) and may correspond to the very late stage 
of the disk evolution in our scenario. Note that this diffusing torus scenario is different from the 
receding torus model in which luminous AGNs may have weak tori due to sublimation (Lawrence 
1991). Generally, the evolution of a torus (or outer disk) will be driven by both diffusion and 
sublimation. 

— Radiatively inefficient accretion: With the decrease of the accretion rate in the second evolutionary 
phase (but before it decreases down to the Bondi accretion rate), the accretion may transit to a 
radiatively inefficient mode (ADAF or ADIOS; Narayan & Yi 1994; Blandford & Begclman 1999). 
The existence of some AGNs whose spectra need to be explained through radiatively inefficient 
accretion modes (but with an accretion rate larger than the Bondi accretion rate) may thus also be 
explained as a consequence of the evolutionary sequence described in this paper. The distribution 
of the normalized accretion rates of nearby low-luminosity AGNs could provide constraints on the 



lr Thc shortness of the characteristic timescale for QSOs can be explained as follows : (i) for a given BH mass, the peak nuclear 
luminosity decreases with increasing characteristic timescale (see eq. 12); (ii) because the local BH mass function decreases 
exponentially at the high-mass end, a long characteristic timescale results in QSO number densities that are much lower than 
the observed ones. Note that the characteristic timescale here is different from the "lifetime" of the nuclear activity phase. The 
value of the QSO lifetime is sensitive to the detailed definition of "nuclear activity" (e.g., the lower threshold of the nuclear 
luminosity for a galaxy to be classified as "active"). The determination of lifetime is more likely to be affected by uncertainties 
of the QSO luminosity function at the faint end. See more discussions in Yu & Lu (2004a). 
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transition from the radiatively efficient accretion mode to the radiativcly in-efficient mode (see 
also Bcgclman & Cclotti 2004). 

— Radio-loud vs. radio-quiet: AGNs/QSOs can be classified into two groups: radio-loud and radio- 
quiet. It has been a long-standing question whether all AGNs pass through a radio- loud phase 
or whether only certain objects develop jets and become radio-loud. Boroson (2002) shows that 
most radio-quiet QSOs have high Eddington ratios, while most radio-loud QSOs have low ratios. 
Since the accretion rate decreases with time as the accretion disk evolves, we propose that from 
the evolutionary point of view, a QSO / AGN appears radio-quiet at its early evolution stage and 
appears radio-loud at its late stage (of course, the detailed appearance is also complicated by 
orientation effects and possibly by BH spins). 

6. Conclusions 

In this paper we have connected the accretion rate distribution of nearby SDSS AGNs (z < 0.3) with the 
accretion physics around massive BHs. The observed distribution of the AGNs (with host galaxy velocity 
dispersion a ~ 70 — 200 km s _1 ) follows a power-law distribution at low accretion rates (mjoin] — 10~ 3 — 0.1), 
which is consistent with a scenario in which their mass accretion rates are declining with time in a power-law 
form [M oc t~ 7 with 7 = 1.26 ± 0.1) and the accretion process follows a self-similar evolutionary pattern, 
as simple theoretical models predict. The transition timescale of the accretion to the self-similar phase 
(estimated from the turnover of the observational power-law distribution at high accretion rates) is roughly 
consistent with the diffusion timescale of an accretion disk. The results suggest that the observed diversity in 
the Eddington ratios of AGNs (~ 10~ 3 — 1) is a natural result of the long-term evolution of accretion disks. 
Some other issues deserving of further investigation are discussed, such as a full explanation for the accretion 
rate distribution of AGNs with higher velocity dispersions (a > 200 kms -1 ), the long-term evolution of 
accretion disks (including self-gravitating disks), the BH mass versus velocity dispersion relation in AGNs, 
the evolution of BBHs in QSOs/AGNs, the coevolution of galaxies and QSOs/AGNs, and the unification 
picture of AGNs. 
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